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We derive an extension to the quantum regression theorem which facilitates the calculation of 
two-time correlation functions and emission spectra for systems undergoing non-Markovian evolu¬ 
tion. The derivation exploits projection operator techniques, with which we obtain explicit equa¬ 
tions of motion for the correlation functions, making only a second order expansion in the system- 
environment coupling strength, and invoking the Born approximation at a fixed initial time. The 
results are used to investigate a driven semiconductor quantum dot coupled to an acoustic phonon 
bath, where we find the non-Markovian nature of the dynamics has observable signatures in the 
form of phonon sidebands in the resonance fluorescence emission spectrum. Furthermore, we use 
recently developed non-Markovianity measures to demonstrate an associated flow of information 
from the phonon bath back into the quantum dot exciton system. 


I. INTRODUCTION 

Two-time correlation functions are quantities of fre¬ 
quent interest in many areas of physics. This is partic¬ 
ularly true in quantum optics, where correlation func¬ 
tions of the form {A{ti)B{t 2 )) give the field correlation 
properties of an emitting system such as a driven atom, 
and whose Fourier transform gives the measured spec¬ 
trum [T]. If the governing Hamiltonian can be diago- 
nalised exactly, calculation of the two-time correlation 
function is no more challenging than calculating a one¬ 
time expectation value of the form (H(ti)). However, it is 
more often the case that the emitting system is an open 
system, whose dynamics can only be approximated. In 
this case, since the system operators A and B are evalu¬ 
ated at two distinct times, calculation of the correlation 
function given knowledge of system dynamics alone is 
not at first sight straightforward. The quantum regres¬ 
sion theorem, however, gives a prescription of how such 
correlation functions can be related to more readily ob¬ 
tainable system expectation values [2]- A subtle caveat 
of the quantum regression theorem, however, is that it 
applies only to systems undergoing strictly Markovian 
evolution. It requires that the complete density operator 
of the system and environment factorises at all times, 
and that the reduced system density operator obeys a 
time-independent master equation unn]. 

The requirement of Markovian evolution is typically 
fulfilled in the traditional case of atomic quantum op¬ 
tics due to the extremely short correlation time of the 
electromagnetic environment [in ng. However, more 
recent technological advances in the fabrication of ar¬ 
tificial emitters and the engineering of structured envi¬ 
ronments have given rise to systems whose evolution is 
not purely Markovian, yet whose properties are typically 
probed optically. These systems include semiconductor 
quantum dots (QDs), for which Rabi oscillations [T5HT5] . 
resonance fluorescence [151420) . and single photon emis¬ 
sion pTH55] have all been demonstrated. QDs, how¬ 
ever, exist in a solid-state substrate, and interactions 
with phonons and nuclear spins can modify their emission 
properties [numis] and also give rise to non-Markovian 
behaviour |27U3T] . Additionally, for technological appli¬ 
cations, such as indistinguishable and entangled photon 


sources |5M55] . it is often desirable to place artificial 
emitters in structured photonic environments such as in 
photonic crystals or micro-pillar cavities, which also have 
the potential to lead to non-Markovian behaviour. 


Thus, in order to model the optical properties of these 
ever more exotic systems, it is important to establish how 
two-time correlation functions can be calculated for open 
systems undergoing non-Markovian evolution. We note 
that efforts in this direction have been made [SHin], and 
the conditions under which the regression theorem holds 
have been scrutinised [7]. Many of these, however, rely on 
a number of uncontrolled approximations, such as artifi¬ 
cially enforcing time-locality EllS], or assuming a restric¬ 
tive (rotating wave-like) form of the system-environment 
coupling HIS]. Additionally, it is not clear to what ex¬ 
tent non-Markovian behaviour has any measurable opti¬ 
cal consequences in physically relevant systems. 


In this work we use projection operator techniques to 
derive a non-Markovian extension to the quantum re¬ 
gression theorem, valid to second order in the system- 
environment coupling strength, and invoking the Born 
approximation only at a single fixed initial time. The sec¬ 
ond order expansion restricts the theory to weak-system 
environment coupling regimes for which non-Markovian 
behaviour is typically only present for short times, and 
which is usually very challenging to observe. The key 
advantage of the present work, however, is that this 
short-time behaviour is of a two-time correlation func¬ 
tion, whose spectral counterpart corresponds to a con¬ 
crete readily measurable quantity. Specifically, we apply 
our formalism to the relevant case of a driven QD [164424)] , 
and find that the experimentally observed phonon side¬ 
bands in the emission spectra are a direct consequence 
of non-Markovian behaviour, which the standard Marko¬ 
vian treatment fails to capture. Moreover, we confirm 
true non-Markovianity and indivisibility of the underly¬ 
ing dynamical map by demonstrating that the phonon 
sidebands are associated with a flow of information from 
the phonon environment back into the QD system |36j . 
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II. TWO-TIME CORRELATION FUNCTIONS 
AND THE REGRESSION THEOREM 

We begin by introducing two-time correlation func¬ 
tions and the standard (Markovian) regression theorem. 
We consider a system S interacting with an environment 
and wish to calculate two-time correlation functions 
of the form G{t,T) = {A{t + T)B[t)) = Trs_|_£:[A(t -|- 
T)_B(t)x(0)], where A and B are system operators, x(0) 
is the total system-plus-environment state at t = 0, and 
Trs_|_£; denotes a trace over both S and E. For a time 
independent Hamiltonian E[ we have A{t) = W(t)AU{t) 
with U{t) = exp[—liJt] (we set h = 1), and using the 
cyclic property of the trace we find 

G(t,T) = Trs [A A(t,r)], (1) 

where the system operator A(t, r) is given by 

K{t,T)=TrE[U{r)Bx{t)U\r)\, ( 2 ) 

with x(t) = U(t)xiO)U\t). From Eq. § we see that 
calculation of G{t,T) amounts to calculating something 
analogous to the expectation value of A, but with respect 
to the operator A{t,T) rather than the reduced system 
density operator p{t) = TTE[U{t)x{0)U\t)]. For this 
reason we refer to A(t, t) as the reduced effective density 
operator, and p{t) the reduced physical density operator. 

The standard regression theorem proceeds by observ¬ 
ing that the definition of the effective density opera¬ 
tor A{t, t) in Eq. ([^ bares a strong resemblance to 
that of the reduced physical density operator, p{f) = 
Tr£:[[/(t)x(0)C/'^(t)]. As such, if we know the equation 
of motion for the physical density operator with respect 
to t, say dtp{t) = $p(t), then the reduced effective den¬ 
sity operator will obey the same equation of motion but 
with respect to r, and with a modified initial condition, 
namely drA{t,T) = ^A{t,T) and A(t, 0) = Bp{t). We 
will see, however, that this procedure contains a hidden 
assumption that the total physical density operator y(t) 
factorises for all times mail]. 


A. Effective Density Operator Master Equation 
Using Projection Operators 

To see how this assumption arises, and how it can 
be removed, we now derive the quantum regression the¬ 
orem using the projection operator formalism |37H40j . 
This well-established formalism was originally developed 
to calculate physical density operator master equations, 
and our purpose here is to do the same for the effec¬ 
tive density operator, taking particular care to identify 
places where any approximations have different physi¬ 
cal significance. To begin we must establish an inter¬ 
action picture for the total effective density operator, 
which we define as T(t,T) = U{T)Bx{t)U^ {t) , such that 
A(<,t) = Tr£;[T(t,r)]. We write the total Hamiltonian 
= Ho + aEli, where Hq = Hs -t He with Hs and He 
acting exclusively on S and E respectively. We recall that 


the unitary operators U{t) and Uo{t) are defined as the 
solutions to the differential equations idrU{T) = HU{t) 
and idrUo{T) = HqUq{t), and the interaction picture ef¬ 
fective density operator as T(f,r) = Ui{T)Bx{t)U\{T) 
with Ui{t) = Uq{t)U{t). From these definitions we find 

drf{t,T) = -ia[H/(r),T(f,r)] = aC{T)f{t,T), (3) 

where Hi{t) = Uq{t)HiUo{t) and the Liouvillian £(t) 
is defined to satisfy the second equality. We naturally 
define A(t, r) = Tr£;[T(t, r)], and note that since we 
can write Uq{t) = U s{t)U e{t) with the subscripts in¬ 
dicating whether the operators act on 5 or H we find 
A(t,r) = Us{T)A{t,T)Ug{T). The Schrodinger and in¬ 
teraction picture equations of motion are then related 
through 

drA{t,T) =i[A{t,T),Hs]+Us{T){dr~A{t,T))Ul{T). (4) 

These results demonstrate that the effective density op¬ 
erator has a well-defined interaction picture which facil¬ 
itates the use of the master equation techniques below. 

We now introduce the projection operators V and Q = 
{t — V), which are defined through [38114(1] 

V'f{t,T) = Tr_E[T(t,r)] ® Pr = A(t, r) ® pn, (5) 

where pR is a reference state of the environment. The 
projection operators project the effective density opera¬ 
tor into factorising and non-factorising components, i.e. 
we can write T(t,r) = (7^ -I- Q)T(t,T), where the first 
term factorises by definition, while the second captures 
those components which do not. From these basic def¬ 
initions one can show that = V and Qf = Q, 
while QV = VQ = Q. In what follows we assume 
Et:e[HiPr\ = 0. This is not an approximation, since if 
Eve[Hipr\ = (Hi) ^ 0 we can redefine H'g = Hs + {Hi) 
and H'j = Hi — {Hi) leaving the total Hamiltonian un¬ 
changed, and we then have Tt:e{H'iPe\ = 0 by defini¬ 
tion mma. Provided our reference state is chosen such 
that [He,Pe\ = 0, valid for e.g. thermal states, we find 
Tt:e[H){t)pr\ = 0 which implies VC{t)V — 0. 

Now, our aim is to derive an equation of motion 
for the factorising part of the effective density operator 
VT{t,T), from which we can readily obtain A(t,r) = 
Tr£;[7^T(t, r)], and using Eq. Q calculate the two-time 
correlation function. Eollowing Ref. m we act with both 
V and Q on Eq. (§ yielding two differential equations 
which we must solve simultaneously. Inserting 1 = 7^ -I- Q 
on the right hand side and using VC{t)V = 0 the first of 
these becomes 

drVf{t,T) = aVC{T)Qf{t,T), (6) 

while the second involving drQTitjT) can be formally 
integrated to give 

Qf{t,T) = GF{T,0)Qf{t,0) 

+ a[ dsGE{T,s)QC{s)rf{t,s), (7) 

Jo 
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where Gf(t, s) = exp [a/J ds'Q£(s')] with 
the chronological time ordering operator jST]. To ob¬ 
tain a time-local form, from Eq. ^ we see that we 
can write T(t, s) = Gb{t, s)T{t,T), where Gb{t,s) = 
T_j. exp [—a JJ" ds'£(s')] with T_j. the anti-chronological 
time ordering operator. From Eq. Q we then find 

(1 - E(r))Qt(t,r) = (r, 0)Qt(t, 0) + E{T)rf{t,T), 

( 8 ) 

where E(t) = a/J” dsG_F(T, s)Q£(s)7^GB(r, s). Pro¬ 
vided the inverse of the operator (1 —E(t)) exists, Eq. ^ 
can be solved for QT(<,r). Since we are ultimately 
interested in the weak-coupling limit of the system- 
environment interaction a, and since E(t) contains no 
zeroth order term in a, we assume the existence of such 
an operator, and in solving for QT(t, r) we obtain 

Qt(t,r)=(l-S(r))-iE(r)7^t(t,r) 

+ (l-E(T))-iG;^(T,0)Qt(t,0). (9) 

Inserting this formal solution for the non-factorising com¬ 
ponent of the effective density operator into Eq. ([^ for 
the factorising component we find 


that the physical density operator factorises at all times, 
x{t) « p{t) 0 pn, then QT(t, 0) = 0 and the inhomo¬ 
geneous term vanishes. Analogously, in Eq. ( |14[ ) we see 
that in assuming factorising initial conditions, x(0) « 
p(0) (8) PR^ the inhomogeneous term for the physical den¬ 
sity operator vanishes. In these cases the equations of 
motion for the effective and the physical density opera¬ 
tor become identical, i.e. we have dtPx{t) = JC{t)Vx{t) 
and dr'PT{t,T) = IC{T)VT{t,T). We conclude that we 
must make the Born approximation at all times for the 
standard regression theorem to apply. 

We now turn to the key insight of this work which 
allows us to remove the Born approximation. Since 
i? is a system operator, and assuming [He,Pr] = 0, 
it can be shown that QT(t,0) = BU[t)Qx{t), where 
Uo{t)x{t) = UQ{t)x{t)Ul{t). The object Qx(t) represents 
deviations from factorability of the physical density op¬ 
erator. However, we already have an exact form for this, 
namely Eq. (13). Assuming factorising initial conditions 


only, the second term in Eq. (13) is zero, and using what 
remains in Eq. (10) gives 


drVfit, r) = I\t, T)rm + nr)VT{t, r), (15) 


drVT(t,T) =I{T)Q.f{t,{)) + lC(T)Vf(t,T), (10) 

where we have defined the kernels 

I(t) = aVL{T){t - E(r))-iGF(r,0)Q, (11) 

/C(t) = alP/:(T)(l-E(r))-iE(T)lP. (12) 

These expressions constitute an exact equation of motion 
for the reduced effective density operator, with an inho¬ 
mogeneous term which depends on the physical density 
operator through QT(t, 0) = QBxit). 

For these reasons, it what follows it will be useful to 
also consider the evolution for the factorising and non¬ 
factorising parts of the physical density operator x(0- 
For this purpose we use the projection operator methods 
outlined above, and we find that the derivation proceeds 
in precisely the same manner, the only difference being 
that the time argument r is replaced with t and the initial 
condition is {V -I- Q)x(O)- In exact analogy with Eq. ([^, 
we find that the non-factorising part has solution 

+ (l-E(<))-iGf(t,0)Qx(0), (13) 

leading to the equation of motion 


where the new inhomogeneous term is gi ven by T' {t, r) = 
T{T)QBUQ{t){t — E(t))“^E(t)7^. Eq. (15) is an exact 


equation of motion for the reduced effective density op¬ 
erator, in which the inhomogeneous term depends on the 
reduced physical density operator, which obeys the exact 


equation of motion Eq. (14) with Qx(0) = 0- 

Though Eqs. (15) and (14) are exact, calculating ex¬ 
plicit forms for the kernels is difficult. The utility of the 
projection operator approach used here is that it allows 
for a systematic expansion in the system-environment 
coupling strength a. Expanding the kernels appearing in 
Eq. (15) to second order in a and moving back into the 


Schrodinger picture we find 

drl^{t,T)=i[K{t,T),Hs\ +V{K{t,T)) +C{g{t,T)), (16) 
where the effective density operator enters through 

V{A{t,T)) = - f dsTiE[Hi,[Hi{-s),A{t,T)pR]], (17) 
Jo 

and the physical density operator enters through 
C{g{t,T)) = 

rT+t 

- / dsTTE[Hi,B{-T)[Hi{-s),g(t,T)pR]], (18) 


dtvm = mQm +( 14 ) 

with the kernels again given by Eqs. 0 and 


B. Removal of the Born Approximation and the 
Non-Markovian Regression Theorem 


Returning to Eq. (10) for the effective density operator, 


we now consider the inhomogeneous term I(r)QT(t, 0). 
If we were to make the Born approximation, and assume 


with g{t,T) = Us{,T)p{t)ul{T), B{-t) = Us{t)BUI{t), 
and we have absorbed factors of a into the interaction 
Hamiltonians, i.e. aHj —> Hj. Let us review what 
approximations have been made. We assumed factoris¬ 
ing initial conditions, x(0) = p(0) ® pR, and expanded 
the kernels to second order in the system-environment 
coupling strength. From this point onwards no further 
approximations are necessary. Finally, we note that 
p{t) entering Eq. ( |I^ can be found at no additional 
cost since to the same level of approximation we have 
dtp{t) = -i[Hs,p{t)]+V{p{t)). 
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Before proceeding, we note that we can obtain a time- 
independent equation of motion for A(t, r) by making a 
Markovian approximation and let r —> oo in Eq. (161. We 
then hnd C(p(t, r)) = 0 and the inhomogeneous term dis¬ 
appears. In this case the regression theorem is recovered 
since p(t) and A(t,r) obey the same equation of motion. 
Recalling that we also find C{g{t,T)) = 0 when making 
the Born approximation, y(t) « p{t) 0 pn, we conclude 
that in the present context the Markovian approxima¬ 
tion cannot be made without also implicitly making the 
Born approximation. Is the converse also true? Is it 
possible to not make the Markovian approximation by 
leaving the integration limit in Eq. (16) at t, yet at the 
same time make the Born approximation and neglect the 
inhomogeneous term C{g(t,T))7 This is what one would 
obtain naively applying the regression theorem to a non- 
Markovian master equation for the physical density op¬ 
erator. In the following we will see that this approach is 
ill-advised and can give rise to unphysical results. 


III. APPLICATION TO A DRIVEN 
SEMICONDUCTOR QUANTUM DOT COUPLED 
TO ACOUSTIC PHONONS 

We now use our results and consider a driven semi¬ 
conductor QD in a non-Markovian acoustic phonon en¬ 
vironment [HIllIlT]. The QD is described by ground 
and single exciton states \g) and |e), and the laser by a 
constant Rabi frequency Q and detuning <5. In a rotating 
frame, and within the dipole and rotating wave approxi¬ 
mations the Hamiltonian is given hy H = Hs + Hj + He, 
with Hs = 5a^a+{VL/2){a^-\-a), Hj = 9 k{b\+bk) 

and He = J2k'^kbl,bk, where a = \g){e\, a phonon with 
wave vector k and frequency Wfc is described by cre¬ 
ation and annihilation operators 6^ and b^, and we 
take a thermal state for the phonon environment pn = 
exp[—HE/kBT]/TT[exp[—HE/kBT]], with T the sample 
temperature. The exciton-phonon interaction is charac¬ 
terised by coupling constants g^, which ultimately enter 
only through the spectral density J(w) = 9k^i^~^k)- 

Eor coupling to longitudinal acoustic phonons we can 
take the form J{uj) = 77 exp[—(w/wc)^], with rj the 
QD-phonon coupling strength, and Wc the cut-off fre¬ 
quency, whose inverse gives the memory time of the en¬ 
vironment [27]. We tune the laser to the phonon-shifted 
QD transition frequency, 5 = dujJ{uj)/u!, set Q = 
0.12 ps“^, and use the realistic parameters 77 = 0.03 ps^ 
and Wc = 2.2 ps“^, with T = 4 K. The steady state first 
order correlation function of the QD emission is = 

limt_,,oo(cr^(t-|-T)cr(<)), which we calculate with Eq. (U^, 
adding a term r(crA(t, r)crl — i{crV, A(t, r)}) with 1/T^ 
100 ps to capture spontaneous emission. Including spon¬ 
taneous emission in this way assumes that the photonic 
environment is strictly Markovian, and is justified fully 
in the Appendix. Having obtained the first order correla¬ 
tion function the incoherent emission spectrum is defined 
as S{Auj) = Re[/o°“ dr( 5 (i)(r) - g(i)(oo))e-*^“^]. 

Fig. 0 shows the real (a) and imaginary part (b) of 
calculated using the Markovian approximation 




Auj (meV) 



t (ps) 


FIG. 1. Real (a) and imaginary (b) parts of the first order 
correlation function, calculated using a Markovian approx¬ 
imation (solid blue), the full non-Markovian theory (dashed 
orange), and the non-Mar kovi an theory but neglecting the in¬ 
homogeneous term in Eq. (18) (green dotted). Plot (c) shows 
the corresponding emission spectrum with the inset showing a 
different scale on which the Mollow triplet can be seen. From 
the main part of (c) it is seen that only the full non-Markovian 
theory correctly captures the phonon sideband at lower en¬ 
ergies. Plot (d) shows the derivative of the trace distance 
between two states evolved from different initial conditions, 
whose positive values for times ~ 1 ps demonstrates backflow 
of information and true non-Markovianity. 


(taking r —)• 00 in Eqs. (17) and (18), solid blue), the 
full non-Markovian theory (dashed orange), and using 
the naive non-Markovian theory (neglecting the inhomo¬ 
geneous term in Eq. (18), dotted green). We see that 
for times less than the environment correlation time of 
^ 1 ps, all three theories predict quite distinct behaviour, 
reflecting the fact that non-Markovian effects are most 
important and these timescales. Plot (c) shows the cor¬ 
responding incoherent emission spectrum, which on the 
inset scale displays the well-known Mollow triplet. From 
the main part of (c), we see that the Markovian theory, 
which predicts no short time oscillations, correspondingly 
predicts no spectral features at large frequencies. The full 
non-Markovian theory, however, predicts a broad side¬ 
band at lower emission frequencies. This sideband is 
well-known experimentally |20l l43ll^ , and is attributed 
to phonon emission, which our theory supports. Thus, 
the phonon sideband in the emission spectrum is a signa¬ 
ture of non-Markovian behaviour. This is a key feature 
of this work; observation of non-Markovian behaviour 
of one-time expectation values typically necessitates ini¬ 
tialising a system in a well-defined state and tracking 
dynamics on very short timescales (ps in this example). 
Steady-state two-time correlation functions, on the other 
hand, capture fluctuations of a system from equilibrium. 
Non-Markovian behaviour of these fluctuations can be 
much more readily observed since their Fourier transform 
corresponds to an emission spectrum [46] . 


We note that while the phonon sideband has been cal¬ 
culated previously, it has so only in the zero driving limit 
Q —>■ 0 where the model becomes exactly solvable and the 
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Mollow triplet is not present miaiiTi. The theory pre¬ 
sented here works for non-zero fl, allowing us to calculate 
the fraction of power emitted into the phonon sideband, 
which for the realistic parameters used here gives ^ 10%, 
in good agreement with recent experiments |20j . 

Interestingly, it can be seen that the naive non- 
Markovian theory predicts a sideband at higher energies, 
in contrast to both intuition and experimental evidence. 
The inhomogeneous term in Eq. (18) which the naive ap¬ 
proach ignores captures deviations of the true state of the 
environment from the reference state pn used in the mas¬ 
ter equation. For the emission spectrum, these deviations 
are important, since we assumed pn to be a thermal state 
with respect to the QD ground state, which is not the cor¬ 
rect initial condition for an emission process. This reveals 
why in neglecting the inhomogeneous term the sideband 
incorrectly appears at higher energies; since it assumes 
the environment to be in equilibrium with respect to the 
QD ground state, it inadvertently gives dynamics which 
correspond more to an absorption spectrum. We note 
that this correspondence is only approximate, and is not 
expected to be a general feature. 

The steady-state correlation function we have calcu¬ 
lated captures fluctuations of the QD about its steady 
state, and our results suggest these fluctuations are non- 
Markovian in nature. In order for this is be confirmed, we 
calculate a non-Markovianity witness in the form of the 
derivative of the trace distance D(p+,/9_) = i|p-|_(t) — 
P-{t)\, where p+(t) and p-{t) are physical density op¬ 
erator states evolved from two different initial states 

P±{0) = ± <^v) with ay = -i |e)(g| -I- i | 5 )(e| [55] ■ We 

are interested here in the evolution of reduced physical 
density operators since these characterise the behaviour 
of physical QD exciton, and as such use the equation 
of motion dtp{t) = —i[HsT p(t)\ +'D{p{t)) (i.e. with¬ 
out inhomogeneous term). A positive derivative of the 
trace distance is interpreted as a flow of information from 
the environment into the system, and is a sufficient con¬ 
dition to prove indivisibility of the underlying dynam¬ 
ical map, both of which can be considered definitions 
of non-Markovianity [TOl |36l EH]- In Fig. [^d) we show 
^D(p+,/ 9_) calculated using the non-Markovian theory 
(dotted, green), and within the Markovian approxima¬ 
tion (solid, blue). We see that our non-Markovian theory 
gives rise to a time interval during which the derivative 
is positive, confirming true non-Markovian behaviour. 


IV. SUMMARY 

We have developed an extension to the quantum re¬ 
gression theorem, valid to second order in the system- 
environment coupling strength, and invoking the Born 
approximation at a single fixed initial time. These re¬ 
sults have been used to demonstrate that phonon side¬ 
bands in the resonance fluorescence emission spectra of a 
QD are a signature of non-Markovian behaviour. In this 
context, it was shown that this non-Markovian behaviour 
is associated with a flow of information from the phonon 
environment back into the QD exictonic system, which is 


a sufficient condition to prove indivisibility of the under¬ 
lying dynamical map. The projection operator method 
used here is an ideal starting point to include higher or¬ 
der system-environment coupling terms, which can in 
some cases lead to an exact resummation [30]. Finally, it 
will be interesting to investigate how the results obtained 
here can be used to optically quantity non-Markovian be¬ 
haviour [TUI 1551 l55H5n| . 


Appendix A: Effective density operator master 
equation for time-dependent interaction 
Hamiltonians 


Here we give an extension to the results provided 
in the main text which facilitates the inclusion of 
time-dependent interaction Hamiltonians. For a time- 
dependent interaction Hamiltonian we can write the com¬ 
plete Schrddiner picture Hamiltonian in the form H (t) = 
Hs + aHi{t) + He- In this case defining an interaction 
picture proceeds analogously as in the main text, and 
the interaction picture equation of motion for the effec¬ 
tive density operator again takes the form of Eq. (§, 
though now we have 

Hi{t) = uI{t)Hi{t)Uo{t), (A1) 


with Hj{t) the Schrodinger picture interaction Hamil¬ 
tonian at time t, and T(f,r) = Ul{T)U{t -I- 
T, f)[Hx(f)][/'^(f-|-r, t)Uo{T) where the time evolution op¬ 
erator satisfies fo) = H{t)U{t, to) with U{to, to) = 

1. For this time-dependent interaction Hamiltonian 
the derivation of the effective density operator master 
equation proceeds precisely as in the main text, and 
we again arrive at the general expression in Eq. (151, 
the only difference being that the implicit occurrences 
of the interaction Hamiltonians are defined through 
Eq. (Al). Expanding to second order in the system- 
environment coupling strength proceeds analogously, 
though some care must be taken when moving back into 
the Schrodinger picture. For a time-dependent Hamil¬ 
tonian the Schrodinger picture equation of motion for 
the effective density again has the form drA{t, t) = 
-i[Hs,A{t,T)] + 'D{A{t,T)) +C{Q{t,T)), though now 


V{A{t,T)) = 

-[ dsTTE[HiiT,0),[Hi{T - s,-s),A{t,T)pii]], (A2) 
Jo 

and the inhomogeneous term is given by 

r _ 

C{p{t)) = -J dsTiE Hi{t,0),B{-t) 

[Hi{t + T- s,-s),Q{t,T)pii] , (A3) 


and we have defined Hi{ti,t 2 ) = Uo{t 2 )Hj{ti)Uo{t 2 )- 
Note that in order to recover the case for a time- 
independent interaction Hamiltonian we simply set the 
first time argument in Hj{ti,t 2 ) to zero. 
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Appendix B: Inclusion of Spontaneous Emission 
within the Markovian Approximation 


Here we give details of how spontaneous emission can 
be included into the effective density operator master 
equation in the context of the quantum dot (QD) ex¬ 
ample in the main text. To so so we consider an op¬ 
tically driven QD coupled to both a phonon and pho¬ 
ton reservoir. Within the dipole and rotating wave ap¬ 
proximations the total Schrodinger picture Hamiltonian 
in a frame rotating at the laser frequency w/ takes the 
form H{t) = Hs + Hu + Hi 2 {t) + Hei + He 2 where 
Hs = 6cr^a + (H/2)(tTl -|- a), Hu = crJ2k 9k{bl + 6fc), 
Hei = J2k^kblbk, while 

Hi2{t) = Y. h,(aate-“'‘ + (HI) 

q 


and He 2 = S5 where parameters with a k sub¬ 

script refer to phonons, while hq is the coupling con¬ 
stant between the quantum dot and photonic mode g, de¬ 
scribed by creation operator aj and frequency Vq. Since 
the total interaction Hamiltonian Hi{t) = Hi i+H j 2 (t) is 
time-dependent we must use Eqs. (A2) and (|A3|), where 


the trace is now taken over both phonon and photon 
degrees of freedom. Assuming that Hn and Hj 2 {t) 
contain no environment operators that act in the same 
Hilbert space (as is the case in our example), one finds 
that provided TteiIHjiPe] = 0, the ‘c ross ’ terms mix¬ 
ing Hji and Hi 2 {t) in Eqs. ( |A2[ ) and (A3) vanish, and 
we can write drA{t,T) = —i[Hs,A{t,T) -|-I?i(A(t, r)) -|- 
Ci{g{t,T)) -I-I?2(A(t,r)) + C 2 {g{t,T)), where Vi and Ci 
cont ain o nly phonon terms, i.e. they are Eqs. ( |A2[ ) 
and (^) with Hj{ti,t 2 ) Ul{ 12 )Hull 0 ( 12 ), and I?2 

and C 2 contain only photon terms, i.e. Eqs. (A2) and 
( A3 ) with Hi{ti,t 2 ) -> Ul{t 2 )Hi 2 {ti)UQ{t 2 )- As in the 
main text we have Uo{t) = Us{t)UE{t) though now 
UE{t) = ex-p[-i{HEi + HE2)t]- 

Let us consider the term in I?2(A(t, r)) in more detail. 
The relevant interaction Hamiltonian can be written 


Hi 2 {t — s,—s) = d{—s)A\—s)e '^^-l-h.c., (B2) 

where d(—s) = q-iHss^^iHss A(—s) = 

Assuming a zero temperature thermal 


state environment for the photons, i.e. pa = priPr 2 
with pri the state of the phonon environment and 
PR 2 = eyi-p[-^Y.qk'qa\aq]/Tv[eyip[-l3Y,qVqa\aq]] with 
/3 -)■ 00 , we find Tr£;[A'l'A'l'(-s)/9ii;] = TrB[AA(-s)/9_R] = 
Tr£;[AtA(—s)pfl] = 0, and we are left with 

T>2(A(t,T)) = - [ dsTr_E[AA'l'(-s)pii;]e*‘^'® 

Jo 

t) — tT(—s)A(t, r)(T^^ -I- h.c. 

(B3) 


We now make a Markovian approximation, with respect 
to the photon environment only, and approximate the 
remaining correlation function as a delta-function, i.e. 
we take Tve[AA^ {—s)pr\ = C,6{s), in which case we find 

X>2(A(t,T)) = r(^crA(t,T)CT'l' - i{crV, A(t,T)}), (B4) 


where T = 2a^C is the spontaneous emission rate. 
Considering now the photonic inhomogeneous term, 
C 2 ig(t,T)), making the same Markovian approxima¬ 
tion for a zero temperature environment results in 
C 2 {g{t,T) = 0 for all times r > 0 of interest owing to 
the integration limits in Eq. (A3). As such, within the 
Markovian approximation for the photonic environment, 
we can simply neglect the photon terms at a Hamiltonian 
level, provided we add a term equal to Eq. (B4) to the 
equation of motion Eq. (16) in the main text. We note 
that approximating the photonic correlation functions as 
delta-functions is expected to be a good approximation 
for quantum dots in free space or in low Q-factor cavities, 
where photon correlation times of ^ 10“^ — 10“^ ps are 
typically orders of magnitude shorter than the phonon 
bath correlation time of ^ 1 ps [Ml . 
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